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Abstract. We study non-Unear rheology of a simple theoretical model developed to mimic 
layered systems such as lamellar structures under shear. In the present work we study a 2- 
dimensional version of the model which exhibits a Kosterlitz-Thouless transition in equilibrium 
at a critical temperature T^. While the system behaves as Newtonain fluid at high temperatures 
T > Tc, it exhibits shear thinning at low temperatures T < Tc. The non- linear rheology in the 
present model is understood as due to motions of edge dislocations and resembles the non-linear 
transport phenomena in superconductors by vortex motions. 



1. Introduction 

Prictional properties of lubricated system are often strongly influenced by non-linear rheology 
of lubricants. Understanding of the physical mechanism of the non-linear rheology is a very 
important basic issue in the science of friction and condensed matter physics in broader scope. 
Examples of the lubricants include various soft matters [l] , glasses and granular systems [21 [3] . 

Non-linear rheology is believed to arise from some combinations of elastic, plastic and viscous 
deformations. However a unified physical understanding of the mechanism is lacking. The 
purpose of the present work is to develop and analyze a simple statistical mechanical model 
which mimics layered systems, such as those with the lamellar structure, under external shear 
stress. We demonstrate that inspite of its simplicity our model exhibits non-trivial rheological 
properties reminiscent of those in real materials. 

The organization of the present paper is as the following. In the next section we define 
our model where we also sketch a close connection between our rheological problem and the 
transport problem in superconductors. In sec. Owe analyze the phase transition in our model 
by a renormalization group theory and Monte Cairo simulations. Then in sec. |3]we analyze the 
flow curve, i. e. relation between the shear-stress and shear-rate in our model using Langevin 
simulations. We analyze the data in terms of a scaling ansatz which is analogous to that for the 
current-voltage relation in superconductors. Finally in sec. [5] we present our conclusions. 
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Figure 1. A schematic picture of a lamellar structure under shear and our lattice model, a) 
Lamellar structures may be formed for example by surfactants immersed in water. By applying 
external shear stress a on the top and bottom planes the lamellar will exhibit elastic, plastic and 
viscous responses, b) Our lattice model mimic such a layered system stacked into the direction 
of the y-axis under external shear cr xy apphed in the direction of the x-axis. The 'phase' 0(j j) 
represents the displacement at lattice site into the x-direction. Notice that the lattice site 
(i, j) represents a fixed position with respect to the laboratory frame. We denote the lattice size 
into the direction of the x-axis as Lx and that of the y-axis as Ly. 



2. Model 

We consider a layered system which is an assembly of flat elastic sheets stacked on top of 
each other as shown in Fig. [TJ Interactions between the elastic sheets consist of two parts 1) 
mechanical coupling which allows both elastic and plastic deformations and 2) viscous coupling 
due to the presence of some solvents such as water. 



2. 1 . Hamiltonian: a simple model with elastic and plastic deformations 

For simplicity we consider a 2-dimensional model in the present work replacing the elastic 
planes by one-dimensional elastic strings. More specifically we consider a 2-dimensional square 
lattice model of size Lx x Ly shown in Fig. [Hb). Extension to a three-dimensional model is 
straightforward. 

We define dimensionless 'phase' variables ^(ij)'s on lattice sites («,j)'s assuming that 
spatial profile of the density is specified as p{i,j) = Po + Pi cos(6'(jj)) with pQ and pi being 
certain constants. The interactions between the phase variables are described by the following 
Hamiltonian, 

H = ^ Xl(^(i+iJ) ~ %>i))^ - -^'Y. cos(%j) - (1) 

Here J and J' denote the strength of the interactions. The 1st term on the r. h. s represents 
the elastic couplings within each elastic layers and the 2nd term is a sinusoidal coupling which 
allows elastic and plastic deformations between adjacent elastic layers. Although the strength 
of the two couplings should be different in general, we choose both of them to be J = J' = 1 
for our convenience. We note that the effective Hamiltonian Eq. ([1]) is a much simplified one 
as compared with more realistic expressions which take into account splay distortions of the 



layer s[l]. 

Apparently the ground state is given by a spatially uniform 9 which corresponds to a 
crystalline structure. Low energy excitations from the ground state are 1) Goldstone modes: 
smooth spacial variation of 9 and 2) plastic deformations due to dislocation-dipoles as shown 
in Fig. [21 An elementary dislocation-dipole consists of a pair of ±1 'charges' on its two ends 
separated by a lattice spacing /q- In a representation of the configuration of the phase variables 
6'(j j) in terms of unit vectors (cos 0(j j), sin6'(j .,■)) as the XY spins, the ±1 charges appear as 
vorticies and anti-vorticies. 




Figure 2. Elementary dislocation-dipoles. If the difference between the two phase variables 
associated with a interaction bond oriented in the y-direction becomes greater in magnitude 
than vr, the bond slips. In an analogy with the electromagnetism, such an an elementary plastic 
deformation can be regarded as creation of a dislocation-dipole of unit length Iq with ±1 changes 
on both ends as shown in a). Note that the dipoles are parallel to the y-axis. The two possible 
orientations of the dipoles correspond to the sign of the phase slips A0. The latters correspond 
to the two types of edge dislocations created in the crystal as shown in b) where the region of 
high densities are marked by lines. 



Periodic boundary condition is imposed in the x-direction while phase variables at the 
top/bottom layers {j = Ly,l) are regarded as to represent 'walls'. We specify their properties 
shortly later. 



2.2. Dynamics 

We model the dynamics by the following Langevin equation for the phase variables ^(jj) and 
their velocity 'y(ij), 



dvc ■\ dH 

j'=j±l 



where m is the effective mass of the phase variables which we choose to be m = 1 for our 
convenience. The 1st term on the r. h. s. of the 2nd equation represents viscous couplings 



between adjacent elastic layers. We choose the bare viscosity to be r/o = 1 for our convenience. 
The 2nd terms is due to the mechanical couplings and the last term is the thermal noise. 

The thermal noise follows the Gaussian distribution with zero mean and the following 

spatio-temporal correlations, 

{^{i,j)it)^{i',j')i't')) = 2r]okBT5{t - t')5i^i'{6jj' - djj'-i - Sjj'+i). (3) 

Note that the correlation between the thermal noise at neighbouring layers is needed due to the 
viscous (dissipative) copling between them. Here fee is the Boltzmann's constant which we put 
to be 1 for our convenience. In the following we use (...) for thermal averages. 

The phase variables 0{ij) on the top (j = Ly) and bottom {j = 1) layers are regarded as to 
belong to 'rigid walls' which are driven into the opposite directions. More precisely we model 
the walls as 0{i,Ly) = and 0(j i) = 9^^ which are expressed via Fourier series. 



% =2^ cos 



(4) 



Here xcm represents the center of mass position of the wall and Iq is the lattice spacing. We 
choose Iq = 1 m. the following. To mimic a 'rough wall' we choose random values for a^s drawn 
from a Gaussian distribution of zero mean and variance 1 while we choose random values for 
(pnS from a uniform distribution between and 27r. Obviously a 'regular wall' can be selected 
as well by choosing a certain values for a„ and (pn- 

In the present work we drive the walls at constant velocities by enforcing, 

^Cm'^'' = ±^^wall(i) (5) 

where it is for I /II. We define the apparent shear-rate 7 as, 

^'wall 



Ly 



(6) 



2.3. Relation to the transport problem in superconductors 

Here let us mention briefly a remarkable connection between our rheological problem and the 
transport problem in superconductors [D]. 

Apparently a very important issue in the problem of superconductivity is the macroscopic 
transport property: how the Ohmic resistance in the normal phase disappear as the 
superconductivity sets-in. Essential macroscopic properties of the superconductivity are 
determined by ordering of the phase of its order parameter. The phase can be identified with 
the phase variable in our model and the effective hamiltonian can be given as ours Eq. ([1]) but 
with the elastic coupling in the x-direction replaced by a sinusoidal one, i.e. the usual isotropic 
XY model. Then our two-dimensional model correspond to a superconducting film [8]. Shortly 
later we discuss similarity and differences between the equilibrium properties of our model and 
the usual XY model. 

A standard model to study macroscopic transport properties in superconductors are the so 
called resistively-shunted-j unction (RSJ) model (see for example [TO]) in which the couping J in 
Eq. ([1]) is regarded as the strength of the Josephson coupling between superconducting grains. 
The mass m and the bare viscosity ryo in our model Eq. ([2]) are regarded as the capacitance 
and the inverse of the so called shunted-resistance between superconducting grains respectively. 
External forces ±1 are applied on the top/bottom layers to mimic in-coming and out-going 
external electric currents. The latters correspond to nothing but the external shear stresses axy 
in our problem. One then measures the voltage drop V induced in the system which corresponds 
to the shear rate 7 in our problem. Thus the I — V (current vs. voltage) characteristic in the 
superconductors corresponds to the axy — 7 (shear-stress vs. shear-rate) relation which is called 
as flow curve in rheology. 



3. Equilibrium Phase transition 

Here let us discuss some essential features of the phase transition in the present model which 
will provide us a useful basis to analyze the rheology of the model. 

First note that our model given by the hamiltonian Eq. ([T]) is similar to the ferromagnetic XY 
spin model. The difference is that the couplings in the x-direction is elastic in our model while 
the couplings are sinusoidal in both x and y directions in the XY model. In the XY models an 
elementary plastic deformation is creation of a pair of ±1 charges, which need not to form the 
specific type of the dipoles shown in Fig. O In superconductors the ±1 charges correspond to 
the quantized vorticies and anti-vorticies. 

It is well known that the XY model exhibits the Kosterlitz-Thouless (KT) transition at a 
finite critical temperature [S]. We now wish to clarify whether our model, which is extremely 
anisotropic, also exhibits a similar phase transition. To this end we set-up a renormalization 
group theory. 

First by taking a continuous limit we obtain an effective model of a scalar field ip{x, y) in the 
2-dimensional space, whose partition function is given as, 



^(VV^)2-/,dCos (^27r(/?J)/o|^) 



(7) 



where (5 = l/ksT is the inverse temperature, Iq is a unit length in the lattice model and //d is 
the fugacity of a dislocation-dipole. Although the above model resembles the the sine-Gordon 
model [6j which can be obtained by taking a continuous limit of the XY model, the argument 
of the cosine function in the 2nd term of the action S is the derivative while it is ip in the 
usual sine-Gordon model. 

We analyze the renormalization group (RG) flow of the temperature T and the fugacity of the 
dislocation-dipoles. The RG flow represented in the T-^u^^ plane is shown in Fig. [3j Remarkably 
it reveals a KT transition [5l E] at a critical temperature Tc as in the usual sine-Gordon model 
in spite of the strong anisotropy in our model. 

In the RG analysis of our model we are forced to follow not only the flow of the temperature 
T and the fugacity fid but also the ratio between the correlation length in the x-direction and 

in the y-direction. We find the ratio behaves non-trivially as, 

In^oc T>T, (8) 



e± V J 

in the high temperature phase. It means that the system tends to order more strongly within each 
elastic layer than between different elastic layers as expected. Note however that ratio converges 
to just a constant as T ^ implying that the anisotropy does not change the universality. 
As in the usual KT transition, the correlation lengths themselves diverge exponentially fast as 

T ^ r+. 
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T > T,. (9) 



where ^4 is a numerical constant. In the whole low temperature phase T < the correlation 
length of the fluctuation remains oo, i.e. the long-range order is absent and the system remains 
critical. 

For later analyses we need to know the precise value of the critical temperature Tc of the 
original lattice model. To find out Tc , we analyzed the relaxation of an auto-correlation function 



C{t) = (1/iV) ^(cos(%,,)(t) - %,,)(0))) 
(hi) 



(10) 




Figure 3. Renormalization group flow of the fugacity /^d of the dislocation-dipoles vs 
temperature T. For the 'bare' temperatures T <Tc, the fugacity is renormahzed to meaning 
that after coarsegraining the entire system can be regarded as an dislocation-free elastic body 
at some effective (renormahzed) temperature. At T > Tc, the flow of //d exhibits an upturn 
meaning that the dislocations are relevant at high temperatures where the elastic layers flow 
smoothly on top of each other. 

starting from a random initial condition at time t = 0. In the cases of 2nd order phase 
transitions, including the KT transition, such an auto-correlation function is expected to decay 

as 

C{t) oc t-^e-^/^ T > Tc, (11) 

in the high temperature phase T > Tc. The relaxation time r is related the correlation length 
^ via 

r ~ To(eAo)" (12) 

where z{> 0) is the dynamical critical exponent and tq is the microscopic time scale associated 
with the microscopic length scale IIq. Since the correlation length ^ diverges in the limit T — > T^, 
T also diverges as well. Then right at the critical temperature Tc, C{t) exhibits a purely power 
law decay t^'^ with some exponent A(> 0). 

In practice we used the heat-bath Monte Cairo (MC) method and simulated relaxations in 
large systems of sizes 320 x 320 and 640 x 640 by which we could observe C (t) without appreciable 
flnite size effects up to t = 10^ MC steps (MCS). By analyzing C{t) at various temperature we 
found Tc/J ~ 1.15 ± 0.03 and A ~ 0.067 ± 0.002. The latter value of A agrees with that found 
in the 2-dimensional XY model [7J suggesting again that the present model belongs to the same 
universality class as the 2-dimensional XY model. 

4. Non-linear rheology 

Now let us discuss the non-linear rheology of the present model based on our results obtained 
by numerical simulations of the Langevin eq. Eq. ([2]) under constant external shear rate 7. 

4.I. Flow curve 

We first examine the flow curve, i.e. the relation between the shear-rate 7 and shear-stress axy 
in the stationary state. To this end we have performed simulations under constant shear-rate 

^ Random initial configurations are realized by choosing S(ij)(0)'s out of a uniform distribution between and 
27r. 
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Figure 4. The flow curve and its scaling plot, a) Shear-stress axy is plotted versus the shear- 
rate 7. We used systems of sizes 40 x 40. By comparing with the data of 60 x 60 we found 
finite size effects are not significant within the range of shear-rates used in the figure, b) The 
scaling plot for the data at T > Tc at shear rates 7 < 0.01 using the scaling ansatz given by 
Eq. ()15p . Here we used d = 2, z = 2 and the correlation length given by Eq. ([9]) with Tc = 1.15 
determined in sec. [3l The only adjustable parameter used in this scaling plot is the numerical 
factor A in Eq. ([9]) which we choose as A = 1.5. 

7 = v^sii/Ly and evaluated the resultant shear-stress axy between adjacent layers, 



In Fig. m a) we show the data of the shear-stress axy measured under various shear-rates 7 and 
at various temperatures T/J in the double logarithmic plot. 

Except for very high shear rate region 7 ^ 0.01 where the contribution of the bare viscous 
coupling parametrized by the rjo becomes dominant, the shape of the flow curve is non-trivial. 

The relation between 7 and axy clearly exhibits the Newtonian fluid behaviour axy = rj^ at 
high temperatures T > Tc under low enough shear rates 7. Apparently the viscosity tj increase 
as the temperature is lowered toward Tc 

At the same time, s/iear-t/iinnm^ behaviour axy oc 7^"" emerges with some exponent a ~ 1/3 
under higher shear rates except in very high shear rate region 7 > 0.01. Remarkably such a 
power-law region extends to lower shear rates as the temperature is lowered. Right at T = Tc it 
appears that the power law behaviour dominates the entire range of the shear rate except the 
very high shear rate region 7 > 0.01. 

^ The effective shear viscosity behaves as rjeff = o-^y/j ~ 'y~°' meaning that the system flows more easily at 
higher shear rates. Such a behaviour is observed in a variety of systems including various soft matters, glasses 
and granular systems. 
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Figure 5. A schematic picture of a collective excitation of dislocation-dipoles. By creating 
elementary dislocation-dipoles coherently in a closed volume, the interior of the volume exhibits 
a uniform phase slip 27r with respect to the rest of the system. 

In the low temperature phase T < Tc, some finite yield stress lim^^o'^ 7^ would emerge 
if some long range order is established. However, long range order is absent in the present 
2-dimensional model as noted before in sec. [3] and the system remains critical in the whole 
temperature range T < Tc. Indeed the power law behaviour a^y oc continues down 

to lower temperatures with a temperature dependent exponent a(T) which increases up to 1 as 

4-2. Scaling law for the flow curve 

The basic features of the flow curve discussed above suggest that there might be a scaling 
law which explains the flow curve in a unified manner. In the following we first rephrase 
a scaling ansatz[9j, which is originally proposed for the non- linear current- voltage relation in 
superconductors, within our context of rheology. Then we examine the validity of the scaling 
ansatz based on our data. In the following we disregard the trivial contribution from the bare 
viscous force which can be neglected closer to Tc under low enough shear-rates. 

Let us consider a cluster of volume ^'^ in which the phase variables are strongly correlated 
with each other. Presumably the life time of such a cluster is of the order of the relaxation time 
~ ''"0(^/^0)^ as given in Eq. (fT2]) . Under a given shear rate 7, the total phase difference across 
the cluster will become AO = jrS^ within the life time r. The cluster will make a phase slip if 
AO ~ 27r, i.e. 

7ro (Pi ~ 27r. (14) 



Jo, 

Such a phase slip event may be viewed as a collective excitation of dislocation-dipoles within 
a volume ^'^ as shown schematically in Fig. O Preseumably such an event takes place as a 
thermally activated process. The energy barrier associated with it will be of order locrxy^'^~^- 
Thus we expect that the probability for the event is a function of lQCJxy£f'~^ /ksT . 
Based on the above observations we propose the following scaling ansatz, 

^ We note however that if we plot the flow curve not in the double logarithmic plot as in Fig.|3]a) but in a linear 
plot (not shown), we would be tempted to conclude that IvaiAf^o a 7^ below Tc because of the very slow decrease 



where the scahng functions /+ and /_ are for the high temperature phase T > and low 
temperature phase T <Tc respectively. 



4.2.1. High temperature phase In the high temperature phase T > Tc, it is natural to expect 
that the system behaves as a Newtonian fluid, correspnding to the Ohmic resistivity, as far as 
the correlation length ^ remains finite. Right at T = Tc, ^ diverges so that dependence on ^ 
must be eliminated. Then a natural scaling ansatz for /+(x) is, 

1+^ (16) 



It means that the macroscopic shear viscosity rj = a/^ exhibits an anomalous scaling close to 

and that the flow curve exhibits a purely power law, shear-thinning behaviour right at T = Tc, 

a oc (18) 

with the shear-thinning exponent, 

2-d + z , , 

At the lower critical dimension d = 2 and upper critical dimension d = 4, the dynamical 
exponent is 2; = 2 [9]. Thus one finds a = 2/3 for d = 2 and a = for 4 respectively. The latter 
implies shear-thinning is absent in d > 4 where mean-field theories hold. 

We have indeed observed a ~ 2/3 at Tc as we noted before (See Fig. [Ha)). Correspondingly, 
in a superconducting film y ~ /'^ behaviour has been observed experimentally 

Now let us examine the scaling ansatz using our data of the shear-stress axy obtained at 
various shear rates 7 and various temperatures T > T^. As shown in Fig. IHb) the scaling ansatz 
explains very well the crossover from the Newtonian fiuid regime to the shear-thinning regime. 

4-2.2. Low temperature phase In the low temperature phase, just below Tc, the system should 
exhibit the shear-thinning again. On the other hand, at lower temperatures the Arrhenius law 
will hold since the plastic deformations take as thermally activated processes. Thus a natural 
scaling ansatz for f-{x) is, 

/-(x)~| (20) 

It means that /inear viscosity, corresponding to the Ohmic resistance, vanishes in the 7 — > limit. 
However note that the yield stress, corresponding to the critical currents in superconductor, 
defined by strictly taking the limit lim^^o o'xy is zero at any finite temperatures due to 
the presence of thermally activated plastic deformations, sometimes called as creep, even in 
crystalline systems. 

In the case of the present 2-dimensional model which exhibit the KT transition, the 
whole temperature range T < Tc is critical in the sense that the correlation length ^ of the 

d-2-z(T) 

fiuctuation remains 00. Thus the critical behaviour f-{x) ~ x t^-i will persist within 
the low temperature phase with some temperature dependent dynamical exponent [10] z{T) 
which decreases down to as T — > 0, i. e. the system exhibits shear-thinning behaviour 
in the whole low temperature phase with the temperature dependent shear-thinning exponent 
a(r) = z{T)/{l + z{T)) (See Eq. 

In our data shown in Fig. [H a), the data at T < Tc indeed exhibit power law behaviour 
(^xy oc (7)^""^^) with temperature dependent exponent 1 — a(T) which decreases as the 
temperature T is lowered. 



5. Conclusions 

In the present paper we analyzed non-linear rheology in a simple theoretical model which 
mimics layered systems such as those with the lamellar structures under shear. More specifically 
we analyzed in detail a 2-dimensional model which we found to exhibit a Kosterlitz-Thouless 
transition at a finite temperature T^- The flow curve exhibits shear-thinning behaviour below T^- 
The flow curve follows very well a scaling ansatz which we obtained by translating the scaling 
ansatz for the non-linear transport in superconductors to that for our rheological problem. 

We wish to report more details of the present work together with some analysis on other 
features of the non-linear rheology in our system, such as an apparent increase of the viscosity 
as the thickness Ly is made smaller than the correlation length ^_|_, shear-banding and stick-slip 
motions, elsewhere. 
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